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Abstract 



We present spectral density reweighting techniques adapted to the analysis 
of a time series of data with a continuous range of allowed values. In a first appli- 
cation we analyze action and Polyakov line data from a Monte Carlo simulation 
on L t L 3 (L t = 2,4) lattices for the SU(3) deconfining phase transition. We cal- 
culate partition function zeros, as well as maxima of the specific heat and of the 
order parameter susceptibility. Details and warnings are given concerning i) au- 
tocorrelations in computer time and ii) a reliable extraction of partition function 
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zeros. The finite size scaling analysis of these data leads to precise results for the 
critical couplings j3 c , for the critical exponent v and for the latent heat As. In 
both cases (L t = 2 and 4), the first order nature of the transition is substantiated. 



2 



1. INTRODUCTION 

Monte Carlo simulations of finite statistical systems at a coupling (3 = (3q generate 
a time series of correlated data such that for appropriate observables / the arithmetic 
average of measurements 

1 N 

KM = N^ fn (L1) 

n 

becomes an estimator of the expectation value < f > (/3q): 

1 N 

f(Po) = </>(#>) = hm -£/„. (1-2) 

n 

Using the spectral density representation of the partition function 

Z{(3) = dS n(S) exp(/3S), (1.3) 

'S'min 

reweighting techniques allow to calculate the estimator f{(3) for all (3 in a sufficiently 
small neighborhood of f3 . Although reweighting techniques have a long history [1-9], 
the extent of practical improvements became only fully realized after the recent work by 
Ferrenberg and Swendsen [7,8]. In this paper we further elaborate on these techniques 
and cast them in a form suitable for practical applications in lattice gauge theories, where 
one has to deal with a time series of continuous variables. In particular, procedures for 
combining ( "patching" ) data from simulations at various f3o values are discussed in some 
detail. Altogether, we have in mind to establish a model analysis which may provide 
useful guidance for the analysis of data from future large scale simulations of the QCD 
deconfining phase transition, like the U.S. teraflop project [11]. Here, we give a reweighting 
analysis for the SU(3) deconfining phase transition. In a subsequent paper we shall present 
our SU(2) results. (From the point of view of reweighting techniques, the SU(2) theory 
turns out to be the more difficult case, because the action density distributions are barely 
distinguishable from Gaussian distributions.) 

Pioneering numerical work on the SU(3) deconfining transition was done some years 
ago [12]. Renewed interest was stimulated by the APE collaboration [13], who raised 
doubts about the first order nature of the transition. This prompted a number of finite 
size scaling (FSS) studies [14-17] with the result that the first order nature of the tran- 
sition was rather unambiguously established. Here, we give the details promised in [17]. 
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Beyond [17] we present an investigation of autocorrelation times and a theoretical analysis 
concerning the numerical calculation of partition function zeros. The latter point enforces 
some corrections of previously stated results. Aside from quantities which are functions of 
the action (energy) [17], we also analyze now the Polyakov line susceptibility. 

The paper is organized as follows: Section 2 introduces the spectral density method as 
used in this paper, section 3 investigates autocorrelations in computer time (for a review 
see [18]) and makes a connection to error bar calculations by binning [19,20], section 4 gives 
our reweighting calculations for the specific heat and a FSS estimate for the latent heat, 
section 5 is devoted to our calculation of partition function zeros and their FSS analysis. 
Beyond earlier approaches [5,21] a consistency argument is developed. In close analogy 
with work in [15] we give in section 6 a reweighting and FSS analysis for the Polyakov loop 
susceptibility. Summary and conclusions are contained in a final section 7. 



2. SPECTRAL DENSITY MC CALCULATIONS 

We consider the SU(3) Wilson action 

S = J2 S p with S p = ^Tr(U p ). (2.1) 
v 

where the sum goes over all plaquettes of the four-dimensional lattice with volume V = 
L t L 3 . The MC simulation provides us with a time series of measurements for the action: 
S n , (n = 0, N) where N is the total number of measurements. We generate our data by 
means of the computer program published in ref. [22], taking measurements (meas.) after 
every sweep. A certain number of initial sweeps are omitted for thermalization (therm.). 
To have a convenient normalization, we monitor the corresponding action per plaquette 
(not to be confused with S p ) 

Sn = ^, (2.2) 

Vp 

where V p = 6V is the total number of plaquettes. Let us first consider a single MC 
simulation at a fixed value (3 = (3 Q . What is the (3 range of validity for a spectral density 
as obtained from the MC simulation at /3 ? Let us assume it is valid for f3 G [/3 m i n , /3 max ]. 
In this range we are able to calculate estimators for quantities of physical interest by 
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N N 

f((3) = - with F = J]/ n exp(A/?S n ), Z = ^exp(A^5 n ) and A0 = (2.3) 

n=l n=l 

For /3 7^ /9 we call this reweighting. For technical reasons (floating point number over- 
flows) one is actually forced to use in practice = Si — S(f3 ) in these formulas. This 
contributes an irrelevant multiplicative factor and makes the arithmetic computationally 
feasible. With f n = exp(i[3 y S n ), the calculation of partition function zeros can be consid- 
ered as a special case (where the numerator may be omitted in the one histogram case). 
For biased quantities the empirical error bar Af of / is conveniently calculated by means of 
the jackknife method and one may also correct for the remaining bias (see [23] for details). 

The relevant interval A(3 = /3 max — /3 m i n (reweighting range) will shrink with increasing 
volume like* 



A/3 «<r(s) ^ 
as 



= const L- p / 2 L~ d/2 = const L~^ v , (2.4) 



where p = ot/v and the last equality made use of the hyperscaling relation [24] p + d = 2/ v. 
We now have to determine the constant in (2.4). Essentially this boils down to the question: 
From which percentage of our total data at (3q do we still expect meaningful results? 
Clearly, this depends somewhat on whether our statistics is large or small, (/-tiles s q of 
our empirical action density with q in the range 0.025 to 0.1 may still be expected to give 
meaningful answers. This has to be converted into a /3-range. Let qi = q and q2 = 1 — q; 
we define (3 min and /3 max by: 

fimin = ftqi and /3 m ax = Pq 2 i 

where j3 q is given by the implicit equation 



KPq) = s 



This equation is solved numerically for f3 q . Figures 1 — 3 illustrate this procedure for the 
4 • 20 3 lattice at f3 = 5.691. Table 1 gives an overview of our SU(3) data including the 
Anin, /?max and s qi , s q2 values for the choice q = 0.025 (up to errors from conservative 
rounding of /? m i n , Anax)- For L t = 4 Figure 4 depicts the A(3(L) ranges versus L. Our 



* Note that cr(s), the standard deviation of s, goes as L p ^ 2 , while |j9=/? ~ L P ■ 
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Table 1: Data and their reweighting ranges 
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subsequent analysis of autocorrelation times shows that for our present data the choice 
q = 0.025 was too optimistic. However, this does not really matter, because the main 
purpose of calculating [/? m i n , Anax] intervals first is to prevent reweighting calculations at 
absurd (3 values. 



In spin systems it is convenient to work with histograms. For lattice gauge theo- 
ries the action varies typically over a continuous range and a histogram method is not 
recommendable for two reasons: 

i) The size of the histogram bin (i.e. of the action interval deemed to constitute a single 
histogram entry) is an extraneous parameter. It is tedious to have to cope with it. 

ii) Whatever the size of the bin, inevitably part of the information contained in the 
original sample gets lost. 



Instead of artificially introducing histograms, it is more convenient to rely directly on the 
empirical time series for the data. This requires to keep all measurements on disk or tape. 
In our present simulations we kept in double precision the spacelike and timelike plaquette 
expectation values and the real and imaginary Polyakov loop values. This amounts to 
up to 4 * 240, 000 Real*8 data per (L, (3q) simulation point, altogether filling up about 
0.2 gigabyte of disk space (in unformatted storage). Consequently, the feasibility of this 
kind of analysis is tightly linked to the recent progress in storage technology and to the 
availability of large disks. 

To cover a larger (3 range one has to patch MC results from runs at different (3q values 
(/?o +1 > Po, i = 1, P)-, whose validity ranges overlap: 

sg 1 > 4 > 4! 1 > 4i- ( 2 - 5 ) 

Various methods can be found in the literature [3,4,8,9]. The two recent discussions [8,9] 
both aim at minimizing the errors in the resultant estimates for n(S). A crucial difference 
is that [9] fixes the needed relative normalizations of the histograms from data in the 
overlap regions only, whereas [8] exploits a self-consistency condition which was previously 
stated in [4] . The approach [8] yields results even if there is no overlap at all, whereas [9] 
cannot be applied in such a situation. For our purposes, results from patches without 
overlap are assumed to be meaningless and applying the self-consistency condition may 
be somewhat dangerous. For histograms with strong overlaps both methods will converge 
towards identical results. 

More generally, it is useful to patch in such a way that the error of the actually cal- 
culated quantity is minimized. This leads to the following very straightforward approach, 
which we present in a form suitable for the time series analysis of our data. The first 
observation is that any combination 

J(/3) = ^i=i ai * wit h weight factors a; = a,(/3) > (2.6) 

is a valid estimator for (/). In the limit of infinite statistics each single histogram would 
yield the correct results. To simplify our considerations we impose the normalization 

p 

Wj = 1 with Wi = ciiZi. (2.7) 

i=i 

This converts equation (2.6) into 
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P F 
/ = with U = ir. (2.8) 

i=l Zl 

This equation makes clear that the optimal choice for the normalized weight factors Wi is 
simply the inverse variance of f i 

1 

and the overall constant is fixed by the normalization condition (2.7). In practical appli- 
cations the exact variances a 2 (/J are unknown and we have to rely on the empirical error 
bars as estimators: 

wi L— . (2.9b) 

(A/,) 2 

Of course, this just means that the standard way to add estimators weighted by their error 
bars is also the most suitable one for combining estimators from MC simulations at various 
(3q values. However, several complications arise which deserve discussion. 

i) Typically, our data exhibit large autocorrelation times (see next section). This limits a 
serious statistical analysis to using twenty jackknife bins. Imposing a 95 % confidence 
limit (more precisely the [0.025,0.975] confidence interval), the x 2 distribution [25] 
implies that (Af i ) 2 /a 2 (f i ) fluctuates in the range [0.58,2.13]. Our experience is that 
the effect of these fluctuations on f{f3) is harmless as long as only data sets from sim- 
ulations at Pq sufficiently close to f3 are included. However, error bar fluctuations may 
become a serious problem if equation (2.9b) is applied blindly. A highly recommend- 
able cut-off is to set W{ = for f3 £ [/? m i n , Pmax]- It m ay be suitable to constrain the 
included data sets even further. For instance by excluding data sets which: a) taken 
alone give unsatisfactory results and b) have (3 located at one of the edges of 
the validity range. 

ii) Once the weight factors are determined, error bars of f(f3) from the combined statis- 
tics are not calculated by the standard rules of error propagation. Instead, new bins 
are formed, each relying on the combined statistics (2.6). If, after binning, autocor- 
relations still existed in the single data sets, they will become further reduced now as 
each new bin combines data from independent simulations. When appropriate, the 
new bins are used to construct jackknife bins in the standard way. 
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iii) For connected estimators f c = f 2 — f , like the specific heat or susceptibilities, it does 
not follow from equation (2.8) that the optimal weight factor is Wi ~ l/a 2 {f c i ). The 

2 2 — 

reason is that one has to calculate / according to / = w ifi) 2 an d not according 

2 2 

to / = ^iWiifi) 2 - But patching f c would require that f 2 and / be calculated 
with identical weight factors. Fortunately however, this problem seems not to be 
too serious either. Weight factors calculated from / i7 f 2 i or f c i should differ little. 
We have already noted that there are substantial statistical weight factor fluctuations 
which, in the reasonable range, do not influence the final result significantly. Therefore, 
we decided in favour of the simplest solution, namely to use the weight factors Wi ~ 
1/(A/J 2 for the calculation of f^((3). 



3. AUTOCORRELATIONS IN COMPUTER TIME 

It has been emphasized in recent literature [18] that one has to control the integrated 
autocorrelation time fi n t, to be sure about the final error bars. However, in a typical 
realistic situation a Monte Carlo simulation may be perfectly reasonable with respect to 
all calculated quantities, including the confidence intervals implied by their error bars. 
And yet, fi n t may remain the only quantity of interest that cannot be calculated reliably 
in such a simulation. It seems that previous investigations did not pay attention to this 
scenario and, therefore, we find it worthwhile to present some details. 

To recall the concepts we first consider results obtained by a Metropolis generation of 
the Gaussian distribution. We generate 131, 072 = 2 17 numbers; due to a finite Metropolis 
stepsize, successive numbers are correlated. The integrated autocorrelation time is defined 
by: 

2r int (n b ) = p(0) + 2f>(i) with p(i) = f S ° - ffi ~ ! )} . (3.1) 

rr( {{so - s){s - s)} 

For rib — > oo one recognizes that 2ri n t is just the ratio of the correct variance <J 2 (s) divided 
by the naive variance (obtained by using all events s n as if they were independent). A 
convenient way to calculate the correct variance is provided by binning [19]. The connection 
between binning and and the integrated autocorrelation time has been emphasized in [20] . 
Let us partition our data set into Nb = 2 17 ~ k bins of length rib = 2 k (Nb = 32 for k = 12) 
and denote the variance of s after the k th binning by crf(s). Note that 

10 





(3.2a) 



whereas 





(3.26) 



Now 




P(0) + 



2(n b - 1) 



2(n b - 2) 



p(2) + ...+ — p(n 6 -l), (3.3) 



and for n& — > oo the approach towards 2T int {nb) follows from the rapid falloff of p(i), i — > oo. 
We introduce the notation (A^s) 2 for the estimators corresponding to cr 2 (s). 

Using the Gaussian Metropolis data, figure 5 compares the increase of the variance un- 
der multiple binning, (Ajts) 2 / ( A s) 2 , with the direct calculation of the integrated autocor- 
relation time. As expected, for sufficiently large n& identical values (ss 3.7) are approached. 
The convergence is somewhat better for the direct calculation of the integrated correlation 
time, whereas the binning procedure is computationally far more efficient. As usual in 
this paper, the error bars of 2T int are calculated with the double jackknife method [23]. 
However, for the error bar of the empirical variance (A^s) 2 we do not use an estimator, but 
assume instead the % 2 distribution with A^ — 1 degrees of freedom for (iV b — l)(Ajts) 2 /(T 2 (s). 
The error bars depicted in the figures are then chosen such that twice their range, seen from 
the center, corresponds to a normal 95.4 % confidence interval. The depicted symbols are 
the centers of these error bars and not the actually calculated estimators (A/ c s) 2 /(Aos) 2 
(due to the asymmetry of the x 2 distribution, these have somewhat lower values - which 
quickly approach the center as Nf, — > oo). If rib is large, the bins do not only become 
statistically independent, but the central limit theorem implies at the same time that they 
become normally distributed (each bin is an average over a large number of original data). 
For the normal (Gaussian) distribution it is well known [25] that the variance is x 2 dis- 
tributed. Therefore, our assumption is justified whenever the statistical analysis is correct, 
i.e. sufficiently large rib values can be reached such that the correlations between the bins 
can indeed be neglected. 

Figure 6 is the analog of figure 5, calculated with our SU(3) action data from the 4-24 3 
lattice at (3q = 5.691. Obviously, satisfactory convergence is not achieved. It is remarkable 
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that direct calculation of the integrated autocorrelation time gives for large n& a much 
noisier estimator than we obtain from multiple binning. The enlargement of the first part 
of figure 6 in its upper left corner demonstrates this most clearly. Presumably binned data 
are favourably (de-) correlated. Consequently, we now rely on multiple binning. For even 
larger rib (small number of bins N b ) also the error bars of o\ ( s) /a 2 ( s) increase rapidly. 
It is a perfectly possible scenario that Nb = 20 bins are already sufficiently independent 
to allow a statistically reliable calculation of the action s (the 2a confidence level of the 
Student distribution with Nb = 20 is almost normal), but a self-consistent calculation of 
the integrated autocorrelation time is impossible. The binning procedure sheds light on the 
reason. Let us assume we have Nb = 20 independent bins. As already mentioned in this 
paper, it follows from the x 2 distribution with Nj, — 1 degrees of freedom that the 95 % 
confidence interval of (As) 2 /a 2 (s) is [0.58,2.13]. In other words, the correct integrated 
autocorrelation time could be two times larger than the estimate or, equally well, only 
half the estimate. (Remember that the error bars given in the figures are half these 
confidence intervals.) To reduce this uncertainty to approximately ±15 %, we need about 
Nb = 400 independent bins. When, like in gauge theories or especially in full QCD, the 
computational effort is large, computer resources may only allow to create about Nb = 20 
independent bins. This may be fully sufficient to calculate quantities of physical interest 
with reliable confidence intervals. However, to compute the integrated autocorrelation 
time reliably requires a factor 20 more computer time. We conclude that the demand to 
calculate the integrated autocorrelation time in each investigation is exaggerated. Instead 
one may have to work under the assumption that the MC statistics is sufficiently large, for 
instance to give about iV 6 = 20 independent bins, and then support this assumption with a 
number of consistency checks. For example, one may perform Student difference tests for 
estimators from independent MC runs (see the subsequent sections). Such an assumption 
of a sufficiently large statistics would be rather similar in spirit to other assumptions (for 
instance about the approach towards the continuum limit) already made for a typical 
lattice gauge theory simulation. 

As one expects, the situation is under better control for our smaller systems. Figure 7 
depicts the results from our 4 ■ 6 3 lattice at /3q = 5.640. While a direct calculation of the 
integrated autocorrelation time remains inconclusive, its estimate from multiple binning 
is possible: 2T int = 100 ± 10 is consistent with all results from k = 9 (N b = 234) on. 
Proceeding from smaller to larger lattices we obtain rough estimates of 2fi nt for all our 
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data sets. These results are included in table 1 by defining the number of independent 
(indep.) measurements as the number of measurements divided by 2fi nt . One has to 
caution that in the case of L = 20 (240,000 original measurements) we can only rely on 
Nb = 14 bins for estimating 2fi nt . This means an error of about ±50 % (around the center 
value) and no consistency checks towards higher k values. The L = 24 estimates may be 
unreliable altogether. However, for our subsequent purposes the achieved overall accuracy 
seems to be sufficient. 



4. SPECIFIC HEAT 

To characterize phase transitions, action (energy) fluctuations are particularly suit- 
able. The variance of the action is related to the specific heat by cy = (3~ 2 a 2 (S)/V . The 
large L finite size behaviour is 

a 2 (s) = (£(s P -4)£(s P -4)) 

\ p p I 

= const L p J2(( S p ~ S p ) 2 ) = const L p V p a 2 (S p ) with p=-, (4.1) 
p 

where a 2 (S p ) is a constant bounded by (5™ m — 5'™ ax ) 2 /4. It is convenient to use the 
action density (2.1): 

a 2 (s) = V~ 2 a 2 (S) = const L p K p -V(S p ). (4.2) 
The exponent p differentiates three interesting classes: 

1. p = for noncritical behaviour, 

2. < p < d for second order critical behaviour, (4-3) 

3. p = d for a first order phase transition. 

Here the given values of p are defined as those obtained in the limit L — > oo. It is convenient 
to introduce the notion of a strong second order phase transition for transitions with p less 
but close to d (for instance p = 2.99 and d = 3). For second and first order transitions, 
the subleading correction to equation (4.3) is assumed to be non-critical. This leads to the 
following FSS fits 

<j 2 (s) = a\L p ~ d + ci2L~ d for a second order phase transition (4.4a) 
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Results from the action. 





• L 3 


A) 


s 


° 2 o(s) 




bias 


/3max 


2 


g3 


5.094 


0.4390 


10) 


3.42 


10" 


-4 


3.44 (08) 10" 


4 


11 % 


5.0917 (15) 


2 


8 3 


5.090 


0.4329 


;i6) 


2.67 


10" 


-4 


2.70 (07) 10" 


4 


35 % 


5.0908 (10) 


2- 


10 3 


5.090 


0.4249 


'26) 


2.05 


10" 


-4 


2.69 (11) 10" 


4 


102 % 


5.0928 (11) 


2- 


12 3 


5.092 


0.4282 


[31) 


2.31 


10" 


-4 


2.56 (07) 10" 


4 


155 % 


5.0928 (07) 


2- 


12 3 


5.095 


0.4429 


'25) 


1.61 


10" 


-4 


2.62 (31) 10" 


4 


62 % 


5.0928 (10) 


4 


4 3 


5.570 


0.52781 


(48) 


2.51 


10" 


-4 


2.698 (66) 10- 


-4 


4% 


5.5446 (32) 


4 


4 3 


5.610 


0.54101 


(34) 


1.84 


10" 


-4 


none 








4 


4 3 


5.640 


0.54822 


(23) 


2.29 


10" 


-4 


none 








4 


6 3 


5.500 


0.49678 


(13) 


5.32 


10" 


-5 


none 








4 


g3 


5.640 


0.53732 


(25) 


5.91 


10" 


-5 


5.91 (12) 10" 


5 


6% 


5.638 (09) 


4 


g3 


5.645 


0.53905 


(20) 


5.99 


10" 


-5 


6.17 (29) 10" 


5 


3% 


5.614 (14) 


4 


g3 


5.660 


0.54322 


(18) 


5.62 


10" 


-5 


none 








4 


6 3 


5.690 


0.55186 


(14) 


4.66 


10" 


-5 


none 








4 


6 3 


5.740 


0.56188 


(12) 


3.47 


10" 


-5 


none 








4 


8 3 


5.600 


0.52471 


(11) 


2.25 


10" 


-5 


none 








4 


8 3 


5.670 


0.54429 


(21) 


2.58 


10" 


-5 


2.59 (07) 10" 


5 


6 % 


5.6713 (36) 


4 


8 3 


5.693 


0.55098 


(23) 


2.35 


10" 


-5 


2.49 (07) 10" 


5 


10 % 


5.6767 (52) 


4 


8 3 


5.720 


0.55787 


(10) 


1.65 


10" 


-5 


none 








4- 


10 3 


5.600 


0.52458 


(07) 


1.11 


10" 


-5 


none 








4- 


10 3 


5.680 


0.54664 


(19) 


1.36 


10" 


-5 


1.406 (50) 10 


-5 


9 % 


5.6876 (25) 


4- 


10 3 


5.693 


0.55119 


(17) 


1.28 


10" 


-5 


1.346 (50) 10 


-5 


6 % 


5.6842 (32) 


4- 


10 3 


5.710 


0.55558 


(14) 


1.02 


10" 


-5 


none 








4- 


12 3 


5.620 


0.52998 


(05) 


6.31 


10" 


-6 


none 








4- 


12 3 


5.681 


0.54603 


(18) 


7.85 


10" 


-6 


9.23 (55) 10" 


6 


23 % 


5.6909 (30) 


4- 


12 3 


5.691 


0.54990 


(27) 


9.35 


10" 


-6 


9.52 (33) 10" 


6 


21 % 


5.6884 (19) 


4- 


12 3 


5.703 


0.55427 


(12) 


6.42 


10" 


-6 


none 








4- 


14 3 


5.682 


0.54637 


(19) 


5.40 


10" 


-6 


6.50 (32) 10" 


6 


28 % 


5.6882 (14) 


4- 


14 3 


5.691 


0.54950 


(30) 


6.99 


10" 


-6 


7.14 (36) 10" 


6 


25 % 


5.6924 (13) 


4- 


14 3 


5.698 


0.55281 


(13) 


4.51 


10" 


-6 


none 








4- 


16 3 


5.683 


0.54598 


(14) 


3.21 


10" 


-6 


4.95 (42) 10" 


6 


44 % 


5.6902 (15) 


4- 


16 3 


5.691 


0.54923 


(28) 


5.16 


10" 


-6 


5.30 (26) 10" 


6 


28 % 


5.6921 (10) 


4- 


16 3 


5.692 


0.54987 


(28) 


5.40 


10" 


-6 


5.47 (26) 10" 


6 


27 % 


5.6916 (10) 


4- 


16 3 


5.697 


0.55245 


(13) 


3.47 


10" 


-6 


none 








4- 


20 3 


5.690 


0.54839 


(22) 


3.32 


10" 


-6 


3.90 (16) 10" 


6 


54 % 


5.6918 (6) 


4- 


20 3 


5.691 


0.54912 


(28) 


3.88 


10" 


-6 


4.01 (15) 10" 


6 


42 % 


5.6915 (6) 


4- 


20 3 


5.692 


0.54918 


(32) 


3.64 


10^ 


3.90 (24) 10" 


6 


43 % 


5.6929 (7) 


4- 


24 3 


5.691 


0.54781 


(18) 


1.56 


10" 


-6 


4.4 (1.7) 10" 


6 


46 % 


5.6934 (9) 


4- 


24 3 


5.693 


0.55101 


(18) 


1.93 


10" 


-6 


2.89 (35) 10" 


6 


46 % 


5.6910 (8) 



and 



a 2 (s) = a± + (12L d for a first order phase transition. (4.46) 
The first order fit allows to determine the latent heat by the relation 

As = 2y/ai. (4.5) 

A weak first order transition is a first order transition with a small latent heat. Numerically 
it is difficult to distinguish a weak first order from a strong second order transition. The 
FSS scaling relations (4.4) are realized for f3 = f3 c , where (3 C is the infinite volume critical 
coupling. In practice, the exact j3 c value is normally unknown. But we can construct the 
L dependent series of /3 max (L) values defined by 

cr 2 (s;L)(/3) = maximum for f3 = /3 max (L). (4.6) 

The corresponding a 2 (s, /3 max ; L) is denoted by cr^ ax (s). Eqs. (4.4) are also satisfied for 
this series. In addition, this approach yields a precise estimate of f3 c through the fit 

/W^) = 13c + const L~ d . (4.7) 

Exponentially small corrections to this fit are of relevance for small systems [7], but are 
better neglected [10] if the purpose is to estimate (3 C . This is best done by fitting results 
from sufficiently large systems. 

For notational simplicity we drop henceforth the distinction between estimators and 
exact theoretical definitions. Numerical results for each of our data sets are given in table 2. 
No stable values for a max are obtained when /? is too far from /3 max . Similarly, insufficient 
statistics may cause unstable behavior. The 4 • 24 3 lattice at (3q = 5.691 is presumably an 
example. The discrepancies with [17] are only due to the fact that we now use 20 jackknife 
bins, whereas the error bars in [17] were estimated with respect to 32 jackknife bins (on 
the large systems a too generous number). The results for a"^ ax (s) are calculated by means 
of the double jackknife approach [23] and table 2 also lists the additional correction for 
the bias in % of the statistical error bar. Clearly, statistical error and bias of the entire 
statistics have a similar order of magnitude. For future reference we have also included 
the average action s = s((3q) and the variance ctq(s) at (5q in table 2. 
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Table 3: Standard error propagation for the action results. 





• L 3 


^L(*)> goodness 


Anax, goodness 


# 


2 


12 3 


2.56 (07) 10' 


" 4 , 0.85 


5.0928 (06), 


0.94 


2 


4 


4 3 


2.70 (07) 10" 


-4 _ 


5.5446 (32), 




1 


4 


6 3 


5.95 (11) 10" 


" 5 , 0.41 


5.631 (08), 


0.16 


2 


4 


8 3 


2.54 (05) 10" 


~ 5 , 0.32 


5.6730 (30), 


0.40 


2 


4- 


10 3 


1.38 (04) 10" 


" 5 , 0.40 


5.6863 (20), 


0.41 


2 


4- 


12 3 


9.44 (29) 10" 


" 6 , 0.65 


5.6891 (16), 


0.49 


2 


4- 


14 3 


6.78 (24) 10" 


" 6 , 0.20 


5.6905 (10), 


0.03 


2 


4- 


16 3 


5.32 (17) 10" 


" 6 , 0.57 


5.6916 (07), 


0.57 


3 


4- 


20 3 


3.95 (10) 10" 


" 6 , 0.86 


5.6920 (04), 


0.29 


3 


4- 


24 3 


2.95 (34) 10" 


" 6 , 0.40 


5.6921 (06), 


0.05 


2 



Results are calculated from table 2 according to standard error progagation. The last 
column shows how many data sets were combined. 

Using standard error propagation, we combine in table 3 for each lattice size the 
a max( s ) an d Anax estimates of table 2. Besides cr^ ax (s) and /3 max the goodness of fit [26] 
is also listed. It is defined as the likelihood that the discrepancies between the estimates 
of table 2 (lattice size fixed) is due to chance. In case of two data sets we use the standard 
Student test [25] with N = 20 data. For more than two data sets we rely on x 2 and 
assume a Gaussian distribution. If the assumptions are correct (of course there are slight 
corrections), the goodness of fit is a uniformly distributed random number between zero 
and one. Except for the goodness of fit for /3 max from the 4 • 14 3 and 4 • 24 3 systems, all 
other values are reasonable. We tend to attribute the bad fit for the 4 • 14 3 to statistical 
fluctuations (after all we have a rather large number of systems), whereas a closer inspection 
(see below) of the 4 • 24 3 data sets gives rise to the suspicion that the assembled statistics 
is insufficient in this case. 

Table 4 combines the action results by means of the patching method outlined in 
section 2. Within their statistical errors, the corresponding estimates of tables 3 and 4 are 
fully compatible. Only for the 4 • 24 3 lattices does patching reduce the error significantly. 
Figure 8 depicts the cr 2 (s,/3) reweighting calculation for the single 4 • 24 3 lattices, whereas 
figure 9 shows the patched result. From these figures it is clear that the real improvement 
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Table 4: Patching for the action results. 





• L 3 


Cx(s) 


/3max 




#: ' 


weights. 




2 • 


12 3 


2 


55 


(07) 

V / 


10" 


-4 


5 


0928 


(07) 

V / 


2: 


0.67, 


0.33. 






4- 


4 3 


2 


71 


(07) 


10' 


-4 


5 


5439 


(34) 


2: 


0.76, 


0.24. 






4- 


g3 


6 


03 


(13) 


10' 


-5 


5 


.624 


(11) 


2: 


0.36, 


0.64. 






4- 


6 3 


6 


03 


(16) 


10' 


-5 


5 


.620 


(11) 


3: 


0.30, 


0.53, 





17. 


4- 


8 3 


2 


54 


(06) 


10' 


-5 


5 


6719 


(30) 


2: 


0.61, 


0.39. 






4- 


10 3 


1 


37 


(03) 


10' 


-5 


5 


6852 


(18) 


2: 


0.45, 


0.55. 






4- 


12 3 


9 


38 


(29) 


10' 


-6 


5 


6896 


(15) 


2: 


0.34, 


0.66. 






4- 


1 4 3 


6 


70 


(25) 


10' 


-6 


5 


6901 


(09) 


2: 


0.48, 


0.52. 






4- 


14 3 


6 


68 


(20) 


10' 


-6 


5 


6895 


(07) 


3: 


0.39, 


0.44, 





17. 


4- 


16 3 


5 


24 


(15) 


10' 


-6 


5 


6916 


(07) 


3: 


0.17, 


0.41, 





43. 


4- 


16 3 


5 


26 


(14) 


10' 


-6 


5 


6913 


(06) 


4: 


0.15, 


0.35, 





36, 0.14. 


4- 


20 3 


3 


93 


(09) 


10' 


-6 


5 


6920 


(04) 


3: 


0.38, 


0.35, 





27. 


4- 


24 3 


3 


35 


(17) 


10' 


-6 


5 


6919 


(13) 


2: 


0.42, 


0.58. 







The last column gives the number of patched data sets and the relative weights, ordered 
by increasing (3q. When the numbers match, the combined data sets are identical with 
those processed in table 3. Otherwise, the validity ranges of table 1 make it clear which 
data sets have been added. All error bars are calculated with respect to twenty jackknife 
bins and corrected for the bias. 

Table 5: L t = 4 FSS fits of a^ ax (s) 



L range 


(Standard error 
10% 10 2 a 2 


propagation) 

Q 


(Patching) 
10 6 ai 10 2 a 2 


Q 


4 - 


24 


2.22(09) 


1.26(02) 


1Q -23 


2.22(08) 1.27(2) 


1Q -25 


6 - 


24 


2.39(09) 


1.20(02) 


0.40 


2.42(08) 1.18(2) 


0.25 


8 - 


24 


2.47(10) 


1.17(03) 


0.76 


2.49(09) 1.15(3) 


0.81 


10 - 


24 


2.48(12) 


1.16(04) 


0.64 


2.52(10) 1.14(3) 


0.84 


12 - 


24 


2.39(14) 


1.22(06) 


0.77 


2.46(12) 1.17(5) 


0.91 


14 - 


24 


2.41(18) 


1.21(10) 


0.58 


2.51(14) 1.14(8) 


0.97 



These fits use eq. (4.4b): thus, a first order transition is assumed. 

due to patching is in fact much more impressive than noticeable from the cr^^s) error bar 
reduction in table 4 as compared with table 3 (4 ■ 24 3 lattice) . For our other data sets the 
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Table 6: L t = 4 a 2 (s) FSS fits of /3 max (L) 







(Standard 


error propagation) 


(Patching) 




L range 


Pc 


a 


Q 


Pc 


a 


Q 


4 - 


24 


5.6934(3) 


-09.54(21.0) 


0.20 


5.6933(3) 


-09.60(22.0) 


0.15 


6 - 


24 


5.6934(4) 


-10.00(01.0) 


0.15 


5.6936(4) 


-10.60(01.1) 


0.13 


8 - 


24 


5.6931(4) 


-08.40(01.2) 


0.60 


5.6933(4) 


-09.30(01.2) 


0.58 


10 - 


24 


5.6928(5) 


-06.10(01.8) 


0.98 


5.6929(5) 


-07.60(01.7) 


0.80 


12 - 


24 


5.6927(5) 


-05.60(02.6) 


0.96 


5.6929(6) 


-07.50(02.5) 


0.65 


14 - 


24 


5.6926(6) 


-04.90(03.4) 


0.89 


5.6932(7) 


-09.30(03.1) 


0.66 



The P c are extracted using eqs. (4.6) and (4.7). 

Table 7: L t = 2 FSS fits of a max (s) and /3 max (L) 



L range 


(Standard error propagation) 
10% 10 2 a 2 Q 


(Patching) 
P c a Q 


6 - 12 
8 - 12 
10 - 12 


2.39(07) 2.17(25.0) 0.14 
2.53(10) 0.94(71.0) 0.48 
2.38(23) 3.10(03.1) 


5.0929(07) -0.44(38.0) 0.36 
5.0937(09) -1.43(84.0) 0.61 
5.0928(20) 0.00(03.0) - 



Fits to cr max use eq. (4.4b); the f3 c are extracted using eqs. (4.6) and (4.7). 

reweighting analysis of single runs yields already reasonable a 2 (s,[3) pictures. Therefore, 
only one more illustration of patching: Figure 10 gives the combined o~ 2 (s,P) curve from 
all our four 4 • 16 3 data sets. 

Let us first analyse the L t = 4 data. Table 5 collects FSS fits (4.4b) (which assume 
a first order transition) for the cr max (s) estimates of tables 3 and 4. Obviously, results 
from lattices in the range from L = 8 on are fully compatible with our hypothesis that 
the transition is first order. If we now treat p as a free parameter and perform second 
order fits (4.4a), we get p = 3.08 ± 0.23 for the range L = 8 - 24 (patched). Although 
clearly consistent with p = 3, the accuracy is not very restrictive and the error becomes 
even worse if we restrict the fit to larger lattice sizes. Our best estimate of the latent heat 
is depicted in figure 11 and yields 
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As = (3.16 ±0.06) 10 



-3 



(4.8) 



Here we have used the L = 8 — 24 fit for the patched results, as given in table 5. The 
patched results are preferred on the basis of the arguments given in section 2. Monitoring 
the goodness of fit leads to the choice L = 8 — 24. As soon as the goodness of fit is 
reasonable, there is no statistically significant inconsistency between the smaller and the 
larger lattices. Still, omitting some smaller lattices in such a situation may decrease the 
systematic errors. However, this is not relevant anymore within the statistical accuracy. 
The main effect of omitting the smaller lattices is only an increase of the statistical noise. 

Next, we use the /3 max (£) estimates of tables 3 and 4 as input for the FSS fit (4.7) 
and arrive at the results of table 6. Our best estimate of /3 C , corresponding to (4.8), is 



The line of arguments is similar as for our latent heat estimate. 

The analysis of the L t = 2 data is straightforward. The first order nature of the 
transition is much more pronounced than for L t = 4. For the 2 • 10 3 lattice the time series 
and corresponding action histogram are depicted in figures 12 and 13. This should be 
compared with figures 1 and 2. In both cases we have the scale factor L/L t = 5, but only 
for Lt = 2 is the two-peak structure immediately clear. For the latent heat as well as for 
f3 c the L t = 2 FSS fits are now summarized in table 7. For L = 12 (the only size for which 
we have two data sets), the results of straightforward error propagation and patching are 
identical (see tables 3 and 4). Thus we only need one table to show these fits. Our best 
estimates (from L = 6 — 12) are 



P specific heat 



5.6933 ± 0.0004, (L t = 4). 



(4.9) 



As = (3.09 ±0.05) 10 



-2 



(L t = 2) 



(4.10) 



and 



P specific heat 



5.0929 ±0.0007, (L t = 2). 



(4.11) 
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5. PARTITION FUNCTION ZEROS 



Ref. [27] discusses the physical relevance of partition function zeros. Their numer- 
ical calculation was pioneered in Refs. [5,28,29,21,9,17], of which [21,17] are concerned 
with SU(3) lattice gauge theory. In spin systems the action takes discrete values and the 
partition function becomes a polynomial in exp(/3). Under such circumstances [28,9] the 
Newton-Raphson method is convenient to calculate the partition function zeros. When 
the action takes continuous values a time series analysis is more recommendable and we 
calculate the partition function zeros in two steps: first we scan graphically [29] for the 
separate zeros of the real and imaginary part. Figure 14 illustrates this for our 4 • 20 3 
lattice at f3 = 5.692. Re(Z) = is denoted by the crosses and Im(Z) = by the circles. 
A partition function zero is obtained when the lines cross. Second, to compute the precise 
value for the leading zero, we then iterate with AMOEBA [26] (with starting values in a 
sufficiently small neighborhood of the zero). 

Before we can present our SU(3) results we have to clarify some subtle details. For 
smaller SU(3) lattices we noted in [30] that our empirical action distributions are well 
described by Gaussian fits. However, a Gaussian distribution does not give rise to zeros 
in the complex (3 plane. Nevertheless we have reported zeros in [17]. To resolve this 
paradoxical situation, we first study Gaussian random numbers and proceed then with the 
analysis of our SU(3) action distributions. 

5.1. The Gaussian Distribution 

Assume a lattice gauge theory MC simulation (V p = number of plaquettes) and let 

x = s — s. (5-1) 
For non-critical behaviour the measured probability density for x will be 

/A 1 
-exp(-Ac 2 ) with A = — I = aV p . (5.2) 

By reweighting with (3 = j3 x + i(3 y we obtain the partition function up to a normalization 
factor: 

z(8) [a r +o ° 

z (3) = —yiL = \ — exp(Bx + iCx) exp(-Ax 2 ) alx (5.3a) 



Z(3 ) V 7T 



DC 
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Table 8: Partition function zeros. 





• L 3 


Po 




0° 

Hy 


bias x 


biasy 


A/3™ ax 




2 


g3 


5 


0940 


5.0910 (18) 





03960 


(10) 


3 % 


-6 


% 


0.03150 


0.05050 


2 


8 3 


5 


0900 


5.0905 (11) 





01759 


(36) 


2 % 


-10 


% 


0.01290 


0.02180 


2 • 


10 3 


5 


0900 


5.0927 (11) 





00865 


(14) 


-8 % 


-6 


% 


0.00640 


0.01040 


2 • 


12 3 


5 


0920 


5.0928 (08) 





00502 


(08) 


-4 % 


-9 


% 


0.00080 


0.00500 


2 • 


12 3 


5 


0950 


5.0929 (11) 





00495 


(11) 


23 % 


-28 


% 


0.00350 


0.00570 


4 


4 3 


5 


5700 


5.5500 (03) 





09800 


(04) 


-2 % 


-5 


% 


0.05000 


0.10800 


4 


4 3 


5 


6100 


5.5600 (05) 





14900 


(07) 


-17% 


2 


% 


0.0000* 


0.1230* 


4 


4 3 


5 


6400 


5.6070 (06) 





18500 


(06) 


-7 % 


7 


% 


0.0000* 


0.1170* 


4 


6 3 


5 


5000 


none 


















4 


g3 


5 


6400 


5.6540 (10) 





05200 


(25) 


104 % 


-140 


% 


0.04200 


0.06500 


4 


g3 


5 


6450 


5.6560 (05) 





07570 


(64) 


-50 % 


-22 


% 


0.0000* 


0.0660* 


4 


g3 


5 


6600 


5.6420 (07) 





07840 


(47) 


10 % 


35 


% 


0.0000* 


0.0670* 


4 


6 3 


5 


6900 


5.6450 (08) 





08010 


(80) 


-12 % 


-11 


% 


0.0000* 


0.0620* 


4 


g3 


5 


7400 


none 


















4 


8 3 


5 


6000 


none 


















4 


8 3 


5 


6700 


5.6747 (23) 





04660 


(27) 


-2 % 


-12 


% 


0.0000* 


0.0410* 


4 


8 3 


5 


6930 


5.6791 (33) 





04980 


(42) 


5 % 


-16 


% 


0.0000* 


0.0420* 


4 


•8 3 


5 


7200 


none 


















4- 


10 3 


5 


6000 


none 


















4- 


10 3 


5 


6800 


5.6889 (14) 





03010 


(18) 


8 % 


-11 


% 


0.0000* 


0.0270* 


4- 


10 3 


5 


6930 


5.6864 (55) 





02800 


(81) 


7 % 


-7 


% 


0.0030* 


0.0270* 


4- 


10 3 


5 


7100 


none 


















4- 


12 3 


5 


6200 


none 


















4- 


12 3 


5 


6810 


none 


















4- 


12 3 


5 


6910 


5.6896 (17) 





02030 


(07) 


6 % 


-3 


% 


0.0000* 


0.0180* 


4- 


12 3 


5 


7030 


none 


















4- 


14 3 


5 


6820 


5.6886 (18) 





01430 


(09) 


-2 % 


-13 


% 


0.0050* 


0.0140* 


4- 


14 3 


5 


6910 


5.6922 (13) 





01380 


(07) 


% 


-13 


% 


0.0000* 


0.0120* 


4- 


14 3 


5 


6980 


5.6859 (23) 





02020 


(28) 


-69 % 


87 


% 


0.0000* 


0.0130* 


4- 


16 3 


5 


6830 


5.6904 (16) 





01010 


(10) 


-11 % 


-15 


% 


0.00810 


0.01060 


4- 


16 3 


5 


6910 


5.6918 (10) 





01010 


(06) 


-2 % 


-13 


% 


0.0000* 


0.0094* 


4- 


16 3 


5 


6920 


5.6917 (10) 





00960 


(05) 


-2 % 


-14 


% 


0.0000* 


0.0092* 


4- 


16 3 


5 


6970 


none 


















4- 


20 3 


5 


6900 


5.6917 (06) 





00554 


(22) 


-5 % 


-16 


% 


0.00230 


0.00570 


4- 


20 3 


5 


6910 


5.6915 (06) 





00527 


(17) 


1 % 


-6 


% 


0.00120 


0.00540 


4- 


20 3 


5 


6920 


5.6929 (07) 





00531 


(28) 


-1 % 


-15 


% 


0.0000* 


0.0051* 


4- 


24 3 


5 


6910 


5.6931 (07) 





00270 


(02) 


-29 % 


-35 


% 


0.00390 


0.00420 


4- 


24 3 


5 


6930 


5.6913 (09) 





00320 


(04) 


32 % 


-49 


% 


0.00280 


0.00390 



with (defining b x and b y ) 



B = V p ((3 x -(3 )=:V p b x and C = V p (3 y =: V p b y . 
Integration over x gives 



(5.36) 



z{j3) = exp 



B 2 -C 2 



exp i 



BC 



AA ) ~" r \ 2A 
We have zeros of the imaginary part for 



cos 



BC\ 



2A 



+ i sin 



BC\ 



2A 



and of the real part for 



2r77T A 

B = and C= — (n = l,2, 

B 



C= ^M, („ = 0,1,2,...). 



Rewriting equations (5.5-5.6) in terms of 6 X , 6 y we obtain zeros of Im(Z) for 



6 X = (6 y arbitrary) and b y 
and of Re(Z) for 



2nnA 2mva 



(V p ) 2 b 2 



V p b x 



K=^ + ^ A = (2n l^ a , (. = 0,1,2,..). 



(^) 2 & a 



K6 X 



The variance a 2 (\z(/3)\) is easily calculated to be 



B< 



a 2 {\z{(3)\) = exp(^-)-K/3)| 2 . 



(5.4) 



(5.5) 



(5.6) 



, (n=l,2,..) (5.5') 



(5.6') 



(5.7) 



Assume that a numerical calculation has generated N independent data with the proba- 
bility density (5.2). We trust with an about 84% (one sided!) confidence level that the 
calculation of z{j3) will not produce artificial zeros in the (B, C)-range determined by 



(5.8) 



Defining 



X 



exp 



AA 



and Y = exp 



-C 2 
AA 
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the equality of (5.8) becomes 



X = (N + l) 1/2 Y, where (iV + l)" 1/2 < Y < 1. (5.9) 

The argument is only approximate, since numerical results within this confidence radius 
may have some degree of independence, which is difficult to assess. Here they are just 
treated as one event. 

To give a numerical example, we take A = l/(2a 2 ) = 80,000 and V p = 6 • 4 • 14 3 = 
65, 856. We use a Gaussian pseudo random number generator and generate MC data 
according to the probability density (5.2). Subsequently, a reweighting analysis is done 
to determine the zeros. For 1,200 and 76,800 independent data, respectively, figures 15 
and 16 compare the exact values for the leading zeros of Im(Z) and Re(Z) with those 
from the simulation. Using equation (5.9) the expected range of validity for the numerical 
results is also indicated and found to be respected by the data. Namely, the apparent 
crossings (zeros of the partition function) are seen to fall outside the confidence radius. 
In the Gaussian case we know for sure that this means they are numerical fakes. For our 
SU(3) data, we shall therefore have to reject any zeros which fall outside the confidence 
radius. 

Table 9: Patching of Partition Function Zeros. 



L t 


• L 3 




/? 


[/r n ,/3r x ] 


^max 


# : weights. 


2 • 


12 3 


5.0928 


0.00501 (7) 


[5.090,5.095] 


0.0060 


2: 0.63, 0.37. 


4 


4 3 


5.5520 (3) 


0.12300 (6) 


none 


0.117* 


3: 0.74, 0.24, 0.02. 


4 


g3 


5.6500 (7) 


0.07800 (5) 


none 


0.073* 


4: 0.29, 0.37, 0.32, 0.02. 


4 


8 3 


5.6740 (2) 


0.04500 (2) 


none 


0.042* 


2: 0.68, 0.32. 


4- 


10 3 


unstable 


unstable 


none 


none 


2: 0.52, 0.48. 


4- 


14 3 


5.6890 (2) 


0.01450 (6) 


[5.687,5.691] 


0.0156 


3: 0.26, 0.44, 0.30. 


4- 


16 3 


5.6915 (7) 


0.00990 (4) 


[5.688,5.693] 


0.0106 


3: 0.09, 0.45, 0.46. 


4- 


20 3 


5.6921 (5) 


0.00550 (2) 


[5.690,5.694] 


0.0058 


3: 0.27, 0.42, 0.31. 


4- 


24 3 


5.6928 (6) 


0.00290 (2) 


[5.689,5.695] 


0.0044 


2: 0.45, 0.55. 



We patch those data sets for which results are also reported in table 8. The weights are 
in order of increasing f3 . Instead of A/?™ ax we report now [/5™ in , /3™ ax ] as f3 is no longer 
unique. 
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5.2. SU(3) results 

For single runs our leading partition function zeros are collected in table 8. To estimate 
whether they are inside or outside the confidence radii defined by equation (5.9), we use 
the estimated number of independent measurements from table 1 (for instance 44 for the 
4 • 24 3 lattices) and (Tq(s) from table 2 as width of a fictitious Gaussian distribution. This 
leads to the A/3™ ax and /3™ ax values reported in table 8. An example of a zero and its 
confidence radius is given in figure 17 (also showing the definition of A/3™ ax and /3™ ax ). 
An asterix in table 8 indicates that the estimated zero lies actually outside the radius 
of confidence. We see that for the 4 • L 3 lattices most results have problems. The issue 
is actually quite subtle as repetitions with similar statistics lead to reproducible results. 
The reason is that for f3 x = /3q and a Gaussian distribution, Z(f3 y ) falls off exponentially 
with increasing f3 y . As soon as the statistical noise becomes large enough this leads with 
certainty to a fake crossover of real and imaginary zeros, as illustrated in figures 14 and 16. 
Upon a closer inspection of table 8 one may argue that the 2 ■ 12 3 , [3q = 5.092 data set has 
also a problem. However, for the L t = 2 distributions we have a pronounced double peak 
structure and the use of cr (s) from table 2 is not justified. Each of the single Gaussians 
has a much smaller width, leading to confidence radii larger than those reported in table 8. 

To rescue some of our estimates for the 4 • L 3 lattices, we appeal to our patching 
method. The confidence radii for the patched results are estimated by iterating the equa- 
tion 



The final results are collected in table 9. Whereas for lattices of size 4 • 4 3 to 4 • 12 3 we still 
have no conclusive results, we can now determine the leading zero on lattices 4 • L 3 with 
L > 14. For these lattices the FSS fit 



p 



p 




(5.10) 



L x l v with u = e 



(5.11) 



gives 



v = 0.35 ±0.02, (L t = 4) 
24 



(5.12) 



with a goodness of fit Q = 0.26. This is consistent with v = 1/d, i.e. with a first order 
transition. Fitting /3^.(L) with the FFS formula (4.7) yields another estimate of the infinite 
volume critical point: 

^zeros = 5.6934 ± 0.0007, (L t = 4). (5.13) 

The fitted range is L = 14 — 24 and the goodness of fit is Q = 0.73. 

Compared with the L t — 4 situation, the L t = 2 case is unproblematic. All our zeros 
from L = 6 to L = 12 are within the radii of convergence and allow the FSS fit (5.11) 

v = 0.332 ±0.004, (L t = 2) (5.14) 

with the acceptable goodness of fit Q = 0.12. Fitting only the range L = 8 — 12 gives 
v = 0.323(6) with the goodness of fit Q = 0.30. The result (5.14) confirms (with less 
effort) the pioneering work of Karliner et al. [21], who reported v = 0.331(6) from their 
more complicated constrained MC calculation of partition function zeros. Our L = 6 — 12 
estimate of the critical is 



^zeros = 5.0930 ± 0.0007, (L t = 2). (5.15) 

with a goodness of fit Q = 0.42. We see that the /3 c zeros estimates are well consistent with 
the /5 c s P ecific heat results of section 4. 

Table 11: Polyakov Susceptibility by standard error propoagation. 





• L 3 


L~ 


3 Xr 


nax(P), 


goodness 


/3max 5 


good 


ness 


# 


2- 


12 s 


1 


95 


(05) 10" 


-l 


0.83 


5.0926 


(06), 


0.67 


2 


4- 


6 3 


2 


92 


(05) 10" 


-2 


0.49 


5.6957 


(26), 


0.97 


4 


4- 


8 3 


2 


37 


(06) 10" 


-2 


0.03 


5.6930 


(18), 


0.86 


3 


4- 


10 3 


2 


01 


(06) 10" 


-2 


0.86 


5.6891 


(12), 


0.19 


3 


4- 


12 3 


1 


96 


(09) 10" 


-2 


0.74 


5.6913 


(12), 


0.35 


2 


4- 


14 3 


1 


75 


(09) 10" 


-2 


0.17 


5.6906 


(10), 


0.07 


2 


4- 


16 3 


1 


77 


(07) 10" 


-2 


0.99 


5.6913 


(06), 


0.74 


3 


4- 


20 3 


1 


72 


(06) 10" 


-2 


0.65 


5.6917 


(04), 


0.20 


3 


4- 


24 3 


1 


67 


(31) 10" 


-2 


0.34 


5.6920 


(06), 


0.05 


2 



Results are calculated from table 10 according to standard error progagation. The last 
column shows how many data sets were combined. 
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Table 10: Results from the Polyakov susceptibility. 





• L 3 


A) 


P 


I 


-3 


Xo(P) 




3-w 
A max 


(P) 


bias 




2 


g3 


5 


094 


o 


77480 


(0211) 


1 


44 


10" 


-i 


1.52 (003) 


lO" 1 


17% 


5 


0836 


(25) 


2 


8 3 


5 


090 





62510 


(0378) 


1 


58 


10" 


-1 


1.60 (004) 


10" 1 


38 % 


5 


0894 


(12) 


2 • 


10 3 


5 


090 





38700 


(0700) 


1 


49 


10" 


-1 


1.87 (006) 


10" 1 


122 % 


5 


0936 


(16) 


2 • 


12 3 


5 


092 





84400 


(0680) 


1 


13 


10" 


-1 


1.95 (005) 


10" 1 


177% 


5 


0928 


(07) 


2 • 


12 3 


5 


095 





45420 


(0872) 


1 


79 


10" 


-1 


2.01 (028) 


10" 1 


53 % 


5 


0922 


(12) 


4 


4 3 


5 


570 


o 


27737 


(0296) 


2 


67 


10" 


-2 




none 












4 


4 3 


5 


610 


o 


33716 


(0404) 


3 


34 


10" 


-2 




none 












4 


4 3 


5 


640 





35873 


(0370) 


3 


53 


io- 


-2 




none 












4 


g3 


5 


500 


o 


09542 


(0050) 





26 


10" 


-2 




none 












4 


g3 


5 


640 


o 


21732 


(0437) 


1 


81 


10" 


-2 




none 












4 


6 3 


5 


645 


o 


23119 


(0455) 


1 


99 


10" 


-2 


3.06 (014) 


lO" 2 


8 % 


5 


6950 


(06) 


4 


6 3 


5 


660 





27107 


(0430) 


2 


38 


10" 


-2 


2.90 (007) 


10" 2 


5 % 


5 


6970 


(04) 


4 


6 3 


5 


690 





37507 


(0670) 


2 


87 


10" 


-2 


2.88 


(007) 


10" 2 


6 % 


5 


6960 


(06) 


4 


g3 


5 


740 





47735 


(0653) 


2 


74 


10" 


-2 


3.06 (014) 


lO" 2 


16% 


5 


6940 


(05) 


4 


8 3 


5 


600 


o 


09220 


(0114) 


2 


59 


10" 


-3 




none 












4 


8 3 


5 


670 





22317 


(0830) 


1 


80 


10" 


-2 


2.33 (010) 


lO" 2 


11 % 


5 


6915 


(35) 


4 


8 3 


5 


693 





32431 


(1003) 


2 


50 


10" 


-2 


2.50 (008) 


lO" 2 


12 % 


5 


6937 


(22) 


4 


8 3 


5 


720 





45422 


(0555) 


1 


61 


10" 


-2 


2.12 (012) 


10" 2 


23 % 


5 


6922 


(51) 


4- 


10 3 


5 


600 


o 


06563 


(0066) 


1 


20 


10" 


-3 




none 












4- 


10 3 


5 


680 





21619 


(0834) 


1 


72 


lO - 


-2 


2.047 (093) 10" 2 


12 % 


5 


6901 


(17) 


4- 


10 3 


5 


693 





32961 


(0900) 


1 


86 


10" 


-2 


1.984 (069) 10" 2 


17% 


5 


6867 


(18) 


4- 


10 3 


5 


710 





41699 


(0815) 


1 


33 


10" 


-2 


1.997 (160) 10" 2 


30 % 


5 


6922 


(28) 


4- 


12 3 


5 


620 





05761 


(0052) 





10 


10" 


-2 




none 












4- 


12 3 


5 


681 





15435 


(1002) 


1 


16 


10" 


-2 


2.04 


(025) 


10" 2 


28 % 


5 


6925 


(17) 


4- 


12 3 


5 


691 





27161 


(1581) 


1 


92 


10" 


-2 


1.95 


(009) 


10" 2 


30 % 


5 


6903 


(16) 


4- 


12 3 


5 


703 





40754 


(0718) 





87 


10" 


-2 




none 












4- 


14 3 


5 


682 





15595 


(1196) 


1 


14 


10" 


-2 


1.64 


(011) 


lO" 2 


30 % 


5 


6884 


(15) 


4- 


14 3 


5 


691 





24181 


(1895) 


1 


82 


10" 


-2 


1.87 


(012) 


lO" 2 


31 % 


5 


6920 


(12) 


4- 


14 3 


5 


698 





37534 


(0768) 





73 


10" 


-2 




none 












4- 


16 3 


5 


683 





10277 


(1031) 





68 


io- 


-2 


1.78 


(009) 


lO" 2 


102 % 


5 


6902 


(15) 


4- 


16 3 


5 


691 





22025 


(1992) 


1 


72 


10" 


-2 


1.76 


(012) 


10" 2 


33 % 


5 


6915 


(09) 


4- 


16 3 


5 


692 





25229 


(1937) 


1 


73 


10" 


-2 


1.78 


(013) 


10" 2 


30 % 


5 


6914 


(09) 


4- 


16 3 


5 


697 





36245 


(0793) 





68 


10" 


-2 




none 












4- 


20 3 


5 


690 





16688 


(1729) 


1 


41 


10" 


-2 


1.65 


(009) 


lO" 2 


44 % 


5 


6914 


(05) 


4- 


20 3 


5 


691 





20443 


(2119) 


1 


71 


10" 


-2 


1.76 


(008) 


lO" 2 


55 % 


5 


6913 


(05) 


4- 


20 3 


5 


692 





19617 


(2420) 


1 


62 


10" 


-2 


1.73 


(012) 


lO" 2 


52 % 


5 


6926 


(06) 


4- 


24 3 


5 


691 





08825 


(1466) 





73 


102^ 


2.43 


(085) 


lO" 2 


50 % 


5.6933 


(09) 


4- 


24 3 


5 


693 





31780 


(1379) 





76 


10" 


-2 


1.55 


(033) 


lO" 2 


32 % 


5.6909 (08) 



All error bars are calculated with respect to twenty jackknife bins and corrected for the 



Table 12: Patching for the Polyakov Susceptibility. 





• L 3 


T ~ 3 

X max 


(P) 




/3max 


# Patches: weights 


2 • 


12 3 


1.95 (06) 


10" 


-1 


5.0927 (06) 


2: 


0.68, 0.32. 


4 


g3 


2.95 (04) 


10" 


-2 


5.7007 (26) 


4: 


0.14, 0.36, 0.33, 0.17 


4 


8 3 


2.37 (06) 


10" 


-2 


5.6922 (12) 


3: 


0.37, 0.42, 0.21. 


4- 


10 3 


2.00 (05) 


10" 


-2 


5.6885 (12) 


3: 


0.41, 0.47, 0.12. 


4- 


12 3 


1.95 (08) 


10" 


-2 


5.6910 (13) 


2: 


0.31, 0.69. 


4- 


1 4 3 


1.71 (06) 


10" 


-2 


5.6896 (07) 


3: 


0.40, 0.41, 0.19. 


4- 


16 3 


1.73 (06) 


10" 


-2 


5.6910 (05) 


4: 


0.15, 0.36, 0.36, 0.13. 


4- 


20 3 


1.70 (05) 


10" 


-2 


5.6917 (04) 


3: 


0.39, 0.34, 0.27. 


4- 


24 3 


1.88 (09) 


10" 


-2 


5.6919 (12) 


2: 


0.46, 0.544. 



The last column gives the number of patched data sets and the relative weights, ordered 
by increasing (3q. All error bars are calculated with respect to twenty jackknife bins and 
corrected for the bias. 

Table 13: L t = 4 FSS fits to Xm^(P) 







(Standard error 


propagation) 




(Patching) 




L range 


10V 


a 2 


Q 


10 2 ai 


a 2 


Q 


6 - 


14 


1.78(05) 


2.5(2.0) 


0.17 


1.73(4) 


2.7(2.0) 


0.04 


6 - 


16 


1.76(04) 


2.6(2.0) 


0.21 


1.71(4) 


2.7(2.0) 


0.05 


6 - 


24 


1.74(04) 


2.6(2.0) 


0.33 


1.71(3) 


2.7(2.0) 


0.04 


8 - 


16 


1.68(06) 


3.5(5.0) 


0.73 


1.63(5) 


3.8(5.0) 


0.46 


8 - 


20 


1.68(04) 


3.5(4.0) 


0.87 


1.63(4) 


3.8(4.0) 


0.63 


8 - 


24 


1.68(05) 


3.5(5.0) 


0.94 


1.66(4) 


3.5(4.0) 


0.16 


10 - 


24 


1.68(05) 


3.4(9.0) 


0.87 


1.68(5) 


3.2(8.0) 


0.11 


12 - 


24 


1.65(08) 


4.8(2.3) 


0.84 


1.69(7) 


2.8(1.9) 


0.06 


14- 


24 


1.70(10) 


1.8(4.3) 


0.91 


1.79(8) 


-2.5(3.0) 


0.25 



Eq. (6.2) is used to fit the Xmax(P) in tables 11 resp. 12. 
6. POLYAKOV LINE SUSCEPTIBILITY 

Refs. [15,31] have advocated the FSS analysis of the susceptibility of the projected 
Polyakov line as a good indicator of the order of the SU(3) phase transition. Since we have 
measured and recorded the real and imaginary parts of the Polyakov line along with the 
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Table 14: L t = 4 X (P) FSS fits of /3 max (£) 







(Standard 


error propagation) 


(Patching) 




L range 




a 


Q 


Pc 


a 


Q 


6 - 


24 


5.6914(3) 


0.4(5.0) 


0.20 


5.6907(3) 


1.1(5.0) 


0.005 


8 - 


24 


5.6917(4) 


-0.7(8.0) 


0.31 


5.6910(4) 


-0.2(7.0) 


0.060 


10 - 


24 


5.6921(4) 


-2.8(1.3) 


0.94 


5.6918(4) 


-3.4(1.3) 


0.600 


12 - 


24 


5.6920(5) 


-2.3(2.2) 


0.87 


5.6918(6) 


-3.9(2.3) 


0.440 


14 - 


24 


5.6923(6) 


-4.4(3.3) 


0.98 


5.6923(7) 


-6.7(3.0) 


0.740 



Eq. (4.7) is used to fit the (3 m ax in tables 11 resp. 12. 

Table 15: L t = 2 FSS fits of x max (P) and its /3 max (L) 



L range 


10 2 ai 


a 2 


Q 




a 


Q 


6 - 12 
8 - 12 
10 - 12 


18.8(5.0) 
21.2(8.0) 
20.6(1.7) 


-08.3(01.3) 
-26.0(05.0) 
-19.0(20.0) 


0.0007 
0.7100 


5.0942(08) 
5.0942(10) 
5.0915(26) 


-2.29 (56.0) 
-2.31(97.0) 
2.10(04.1) 


0.53 
0.26 



The susceptibility maxima for L t = 2 are fitted to eq. (6.2) and the corresponding (3 max 
to eq. (4.7). 

plaquette action, for each of our runs, we can apply the procedures discussed in Section 4 
to the spectral-density FSS analysis of this quantity. 

We have a time series of measurements of the lattice average of the Polyakov line 
in the Euclidean time direction, O = V~ x J2 X Q X - These are complex numbers: O = 
Re Q+ ilm O = pe 1 ^ . (In our definition, Q x is larger than in Ref. [15] by the colour factor 
3). 

The projected real part P is computed as 

C 1. P = Re fl for e [-tt/3, tt/3), 

^ 2. P = Re exp(-z27r/3)0 for G [tt/3, tt), (6.1) 
^3. P = Re exp(i27r/3)0 for e [-tt, -tt/3). 

P provides a clear distinction between the Z 3 -symmetric phase of the theory (where 
it is close to zero) and any of the three broken phases (where it is nonzero). Since (p is 
projected out, there is no cancellation due to tunneling between the three broken phases 
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(which makes the time series average vanish on a finite lattice, even above deconfine- 
ment). The susceptibility (variance) of P can now be analyzed exactly like the variance 
of the action was analyzed in Section 4. (To compute the moments of P, we calculate the 
Boltzmann factors from the corresponding action measurements.) The results correspond- 
ing to table 2 are collected in table 10. (Remember that our susceptibilities differ from 
those of Ref. [15] by a factor 9, because of the normalization of fi x .) 

Table 11 shows the results of the standard error propagation analysis, as applied to 
the valid results in table 10. Since the 4 • 4 3 lattices yield no valid results, there is no 4 • 4 3 
entry in table 11. Like for the specific heat, the data sets for L = 14, 24 (Lt = 4) give poor 
estimates of I3 max and the error on L~ 3 Xmax is large for L = 24. Again, patching with 
weights given by eq. (2.9b) (where we put fi = Pi) improves the situation, as can be seen 
in table 12. However, we notice that the patched L~ s Xmax for L = 24 (L t = 4) violates 
the trend in the results for the other L. 

Regarding the FSS behavior of the susceptibility maxima and of the corresponding 
values of the same considerations apply as discussed in Section 4. Assuming a first order 
transition, we fit the data in tables 11 and 12 with the analog of eq. (4.4b): 



Tables 13 and 15 show the results of these fits for L t = 4 and L t = 2 respectively. For 
L t = 4, we see that the patching data (which are of better quality) are more restrictive than 
those obtained by error propagation. Our suspicion about the data L = 24 (presumably 
insufficient statistics) are confirmed by the Q values for the corresponding fits. In addition, 
L = 6 is presumably a too small lattice. The best estimate of the order parameter jump 
A? = should be the one obtained from the patched fit L = 8 — 24 (Q is still 

acceptable, and while it makes sense to ignore the L = 6 result we cannot discard our 
largest lattice L = 24 even though the statistics is poor.) We get 



L- 3 X (P) = ai + a 2 L 



-3 



(6.2) 



AP = (0.26 ±0.04), (L t = 4). 



(6.3) 



For L t = 2 too, the L = 6 data spoil the first order fit; the best estimate is 



AP = (0.92 ±0.18), (L t = 2) 



(6.4) 



from the range L = 8 — 12. 



29 



We can also use the (3 values corresponding to the maxima of the susceptibility, in 
order to obtain (3 C through the fit (4.7). The results of this exercise appear in tables 14 
(for L t = 4) and 15 (for L t = 2). Our best estimates for /3 C SUSC are 



susc 



5.6918 ±0.0004 (L t = 4) 



(6.5) 



(using L = 10 — 24 by patching), and 



susc 



5.0942 ± 0.0008 (L t = 2) 



(6.6) 



using L = 6 — 12. Note that the (3 C fit selects optimal ranges which are different from the 
one we used for the fit to the susceptibility maxima. Obviously, this is allowed in principle 
since the height of the peak and its location in (3 may be independent functions of L. 

While (3 C SUSC for Lt = 2 is seen to be consistent with those estimated from the analysis 
of the specific heat (cf. eq. (4.11)) and of the partition function zeros (cf. eq. (5.15)), 
j3 c snsc for L t = 4 is rather small and becomes only consistent on a two a level (cf. eqs. 
(4.10), (5.13)). This may indicate that with our statistics (presumably due to long time 
correlations) the Polyakov susceptibility is not accurate enough to allow really precise fits. 

7. SUMMARY AND CONCLUSIONS 

Spectral density methods greatly facilitate accurate FSS calculations. One can calcu- 
late pseudocritical couplings precisely and extrapolate towards the infinite volume critical 
coupling. From the specific heat, the analysis of the partition function zeros and the 
Polyakov loop susceptibilities we have obtained three estimates of the infinite volume crit- 
ical (3 which differ somewhat due to remaining systematic errors. In the absence of other 
strong criteria, one may average these estimates, weighted by their error bars and quote as 
the final error the best error bar of the single estimates (one can not use error propagation 
as all results are obtained from the same configurations). In this way we obtain from (4.9), 



(5.13) and (6.3) 



(3 C = 5.6927 ±0.0004, (L t = 4), 



(7.1) 



and from (4.1), (5.15) and (6.4): 



p c = 5.0933 ± 0.0007, (L t = 2), 
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(7.1) 



The achieved accuracy improves estimates from the pioneering literature [12] by one order 
of magnitude in the error bar (translating to two orders of magnitude in computer time). 

Another notable result are the latent heats (4.8) and (4.11), which are consistent with 
independent results by other groups [15,32]. Again the spectral density FSS approach works 
quite well. The possibility to calculate a latent heat and, similarly, an order parameter 
jump self-consistently is in our opinion the strongest argument in favour of the first order 
nature of this phase transition. Whereas for L t = 2 one observes, additionally, a clear 
double peak structure the Lt = 4 transition is fairly weak and a double peak structure 
begins (marginally) only to develop from L = 16 on. It is remarkable that the FSS behavior 
is nevertheless already quite indicative for the first order nature. This seems to be driven 
by the increase of the width of the almost Gaussian looking action density distribution. 

For Lt = 4 the analysis of the partition function zeros has turned out to be more 
subtle than previously anticipated. Nevertheless from L = 14 on the results seem to be 
conclusive and the obtained estimate (5.12) of the critical exponent v is consistent with 
v = 1/d (d = 3) and rounds the picture of a weak first order transition. For L t = 2 
the same analysis is fairly unproblematic and the results of [21] are nicely improved and 
confirmed. The L t = 2 estimate (5.14) for v is, of course, consistent with a first order 
transition. 
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